names(rigal_trends) <- var_temporal_trends
slope <- map_dfr(rigal_trends, ~.x %>%
  select(siteid, linear_slope),
.id = "response"
)
slope_df <- slope %>%
  pivot_wider(names_from = "response", values_from = "linear_slope")

4 Comparison rigal, lm

rigal_trends_df <-  map2_dfr(
  rigal_trends, names(rigal_trends),
  ~.x %>% mutate(variable = .y)
  ) %>%
  select(variable, siteid, linear_slope) %>%
  pivot_wider(names_from = "variable", values_from = "linear_slope") %>%
  mutate(type = "rigal")
slope_comp <- list()
slope_comp$rigal <- rigal_trends_df[order(rigal_trends_df$siteid), ]
slope_comp$simple_lm <- slope_df[order(slope_df$siteid), ]
stopifnot(all(slope_comp$rigal$siteid == slope_comp$simple_lm$siteid))
slope_comp$diff <- tibble(
  siteid = slope_comp$rigal$siteid)
for (i in var_temporal_trends) {
  slope_comp$diff[[i]] <- slope_comp$rigal[[i]] - slope_comp$simple_lm[[i]]
}
old_par <- par()
par(mfrow = c(3, 2))
for (i in var_temporal_trends) {
  plot(
    slope_comp$rigal[[i]] ~
      slope_comp$simple_lm[[i]],
    ylab = "Rigal (polynomial degree 2)",
    xlab = "Simple LM"
  )
  abline(0, 1)
  title(i)
}

par(old_par)
#> Warning in par(old_par): le paramètre graphique "cin" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cra" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "csi" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cxy" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "din" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "page" ne peut être changé
# Compare Rigal and simple lm for few sites
plot_rigal_raw_data <- function(
  site = NULL,
  response = NULL,
  dataset = NULL,
  rigal_trends = NULL
  ) {

  raw_data <- dataset %>%
    filter(siteid == site)
  rigal_resp <- rigal_trends[[response]]
  rigal_resp_site <- rigal_resp[rigal_resp$siteid == site, ]
  plot(raw_data[[response]]~raw_data[["year"]], ylab = response, xlab = "Year")
  abline(rigal_resp_site$intercept, rigal_resp_site$linear_slope, col
  = "red")
  abline(lm(raw_data[[response]]~raw_data[["year"]]), col = "green")
  title(main = site)
}
plot_rigal_raw_data(site = "S2672", response = "log_total_abundance", dataset =
analysis_dataset, rigal_trends = rigal_trends)

slope_comp$diff <- slope_comp$diff %>%
  arrange(desc(abs(log_total_abundance)))
old_par <- par()
par(mfrow = c(3, 2))
for (i in slope_comp$diff[1:6,]$siteid) {
  plot_rigal_raw_data(
    site = i,
    response = "log_total_abundance",
    dataset = analysis_dataset,
    rigal_trends = rigal_trends
  )
}

par(old_par)
#> Warning in par(old_par): le paramètre graphique "cin" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cra" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "csi" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cxy" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "din" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "page" ne peut être changé

5 Polynomial

classification_station <- map2(rigal_trends, names(rigal_trends),
  ~.x %>%
  ggplot(aes(x = shape_class)) +
  geom_bar() +
  labs(x = "Trajectory classification", y = "Nb station",
    title = get_var_replacement()[.y]) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
)


plot_grid(plotlist = classification_station)

p_classif_abun_rich <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[str_detect(names(classification_station),
  "abundance|richness|species")]])
save_plot(
  here("doc", "fig", "classif_trends_abun_rich.png"),
  p_classif_abun_rich
)
p_classif_hill_evenness <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[str_detect(names(classification_station),
  "shannon|simpson|evennes")]])
save_plot(
  here("doc", "fig", "classif_trends_hill_evenness.png"),
  p_classif_hill_evenness
)
p_classif_turnover <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[names(classification_station) %in%
      c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
        "disappearance")]])
save_plot(
  here("doc", "fig", "classif_trends_classif_hill_evenness.png"),
  p_classif_turnover
)
ti <- map_dfr(rigal_trends, ~tabyl_df(x = .x, group = "direction"),
  .id = "response"
)
ti %>%
  filter(direction != "Total") %>%
  select(response, direction, percent) %>%
  pivot_wider(names_from = "direction", values_from = "percent")
#> # A tibble: 18 × 4
#>    response            stable increase decrease
#>    <chr>               <chr>  <chr>    <chr>   
#>  1 total_abundance     75.8%  14.3%    9.9%    
#>  2 log_total_abundance 73.3%  16.5%    10.2%   
#>  3 species_nb          77.7%  13.2%    7.4%    
#>  4 log_species_nb      76.2%  13.2%    6.6%    
#>  5 chao_richness       79.8%  10.6%    8.2%    
#>  6 chao_shannon        76.5%  12.1%    7.7%    
#>  7 chao_simpson        78.4%  10.5%    7.4%    
#>  8 chao_evenness       82.8%  7.5%     5.9%    
#>  9 jaccard             76.0%  3.5%     20.4%   
#> 10 horn                75.9%  4.4%     17.6%   
#> 11 chao                75.7%  4.0%     17.4%   
#> 12 hillebrand          75.3%  4.1%     18.6%   
#> 13 total               73.7%  18.3%    3.0%    
#> 14 appearance          70.2%  15.6%    3.8%    
#> 15 disappearance       64.0%  9.5%     3.0%    
#> 16 evenness            78.8%  9.4%     8.1%    
#> 17 shannon             75.7%  12.2%    8.3%    
#> 18 simpson             77.9%  10.8%    7.7%
#3.3   Log-linear model:logYi=α+βXi+εiIn

#the log-linear model, the literal interpretation of the estimated
#coefficientˆβis that a one-unitincrease inXwill produce an expected increase in
#logYofˆβunits. 
#In terms ofYitself, this
#meansthat  the  expected  value  ofYis  multiplied  byeˆβ.   So  in  terms  of
#effects  of  changes  inXonY(unlogged):•Each 1-unit increase inXmultiplies the
#expected value ofYbyeˆβ.•To compute the effects onYof another change inXthan an
#increase of one unit, call thischangec, we need to includecin the exponent.  The
#effect of ac-unit increase inXis tomultiply the expected value ofYbyecˆβ. So the
#effect for a 5-unit increase inXwould bee5ˆβ.•For small values ofˆβ,
#approximatelyeˆβ≈1+ˆβ. We can use this for the following approxima-tion for a
#quick interpretation of the coefficients:  100·ˆβis the expected percentage
#changeinYfor a unit increase inX.  For instance forˆβ=.06,e.06≈1.06, so a 1-unit
#change inXcorresponds to (approximately) an expected increase inYof 6%
#https://data.library.virginia.edu/interpreting-log-transformations-in-a-linear-model/
hist_linear_plot <- map2(rigal_trends, names(rigal_trends),
  ~.x %>%
  ggplot(aes(x = linear_slope)) +
  geom_histogram() +
  labs(x = "Linear slope", y = "Number of station",
    title = get_var_replacement()[.y]) +
  ylim(c(0, 1000))
)
plot_grid(plotlist = hist_linear_plot)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

type_response <- list(
  abun_rich = str_detect(names(classification_station),
    "abundance|richness|species"),
  hill_evenness = str_detect(names(classification_station),
    "shannon|simpson|evennes"),
  turnover = names(classification_station) %in%
    c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
      "disappearance")
)

map2(type_response, names(type_response),
  function(x, y) {
  p <- plot_grid(
    plotlist = hist_linear_plot[
      names(hist_linear_plot)[x]
      ]
  )
  save_plot(
    here("doc", "fig", paste0("hist_linear_slope_", y, ".png")),
    p
  )

  }
)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> $abun_rich
#> NULL
#> 
#> $hill_evenness
#> NULL
#> 
#> $turnover
#> NULL

5.2 Analysis

5.3 Reproducibility

Reproducibility receipt

## datetime
Sys.time()
#> [1] "2022-01-10 17:38:33 CST"

## repository
if(requireNamespace('git2r', quietly = TRUE)) {
  git2r::repository()
} else {
  c(
    system2("git", args = c("log", "--name-status", "-1"), stdout = TRUE),
    system2("git", args = c("remote", "-v"), stdout = TRUE)
  )
}
#> Local:    main /home/alain/Documents/post-these/isu/RivFishTimeBiodiversityFacets
#> Head:     [34e0a3d] 2022-01-10: Compute avg abundance over first years

## session info
sessionInfo()
#> R version 4.0.5 (2021-03-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 10 (buster)
#> 
#> Matrix products: default
#> BLAS:   /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRblas.so
#> LAPACK: /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRlapack.so
#> 
#> locale:
#>  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
#>  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
#>  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] untb_1.7-4              future_1.21.0           vegan_2.5-7            
#>  [4] lattice_0.20-41         permute_0.9-5           codyn_2.0.5            
#>  [7] janitor_2.1.0           cowplot_1.1.1           rnaturalearthdata_0.1.0
#> [10] rnaturalearth_0.1.0     sf_1.0-4                rmarkdown_2.11         
#> [13] scales_1.1.1            kableExtra_1.3.1        here_1.0.1             
#> [16] lubridate_1.7.9.2       magrittr_2.0.1          forcats_0.5.1          
#> [19] stringr_1.4.0           dplyr_1.0.4             purrr_0.3.4            
#> [22] readr_2.1.1             tidyr_1.1.2             tibble_3.1.6           
#> [25] ggplot2_3.3.3           tidyverse_1.3.0         tarchetypes_0.3.2      
#> [28] targets_0.8.1           conflicted_1.1.0       
#> 
#> loaded via a namespace (and not attached):
#>  [1] colorspace_2.0-0   ellipsis_0.3.2     class_7.3-18       rgdal_1.5-28      
#>  [5] rprojroot_2.0.2    snakecase_0.11.0   fs_1.5.1           rstudioapi_0.13   
#>  [9] farver_2.0.3       listenv_0.8.0      fansi_0.5.0        mathjaxr_1.4-0    
#> [13] xml2_1.3.2         partitions_1.10-4  codetools_0.2-18   splines_4.0.5     
#> [17] cachem_1.0.4       knitr_1.36         polynom_1.4-0      jsonlite_1.7.2    
#> [21] broom_0.7.4        cluster_2.1.1      dbplyr_2.1.0       compiler_4.0.5    
#> [25] httr_1.4.2         backports_1.2.1    assertthat_0.2.1   Matrix_1.3-2      
#> [29] fastmap_1.1.0      cli_3.1.0          htmltools_0.5.1.1  tools_4.0.5       
#> [33] gmp_0.6-2.1        igraph_1.2.6       gtable_0.3.0       glue_1.5.1        
#> [37] Rcpp_1.0.6         cellranger_1.1.0   jquerylib_0.1.3    vctrs_0.3.8       
#> [41] nlme_3.1-152       xfun_0.28          globals_0.14.0     ps_1.6.0          
#> [45] rvest_0.3.6        lifecycle_1.0.1    MASS_7.3-53.1      Brobdingnag_1.2-6 
#> [49] hms_1.1.1          parallel_4.0.5     yaml_2.2.1         memoise_2.0.0     
#> [53] sass_0.3.1         stringi_1.7.6      highr_0.9          corrplot_0.84     
#> [57] e1071_1.7-4        rlang_0.4.12       pkgconfig_2.0.3    evaluate_0.14     
#> [61] labeling_0.4.2     processx_3.5.2     tidyselect_1.1.1   parallelly_1.23.0 
#> [65] bookdown_0.24      R6_2.5.1           generics_0.1.0     DBI_1.1.1         
#> [69] pillar_1.6.4       haven_2.3.1        withr_2.4.3        mgcv_1.8-34       
#> [73] units_0.6-7        sp_1.4-5           modelr_0.1.8       crayon_1.4.2      
#> [77] KernSmooth_2.23-18 utf8_1.2.2         tzdb_0.2.0         grid_4.0.5        
#> [81] readxl_1.3.1       data.table_1.13.6  git2r_0.29.0       callr_3.7.0       
#> [85] reprex_1.0.0       digest_0.6.27      classInt_0.4-3     webshot_0.5.2     
#> [89] munsell_0.5.0      viridisLite_0.3.0  bslib_0.2.4